
addpath ../Generator
addpath ../Generator/GrTheory
addpath ../PDCO
addpath ../CRNTSolvers

mass = [0.0003; 0.5002; 0.5002; 0.4998; 0.4998];
%Generate a random single linkage class and try 
%to converge to the fixed point of the mapping
if(~isvarname('Y'))
		m = 9;
		n = 15;
		r = 3;
		l = 2;
		%Set the seed to be able to reproduce the experiment
		RandStream.setDefaultStream(RandStream('mt19937ar','seed',0));
		%Generate a random set of n complexes on n species, each with at most 3 species
		Y = YGenerator(m,n,r);
		%Generate the strongly connected graph with a single linkage class`
		Ak = AkGenerator(n,0.2,l);
        mass = ones(m,1);
else
		m = size(Y,1);
		n = size(Y,2);
end
A = Ak' - diag(diag(Ak));
D = -diag(diag(Ak));

YAt = Y*A';
YD = Y*D;
W=[];
for i = 1:n
    for j = 1:n
        if A(i,j)~=0;
            W = [W Y(:,i)-Y(:,j)];
        end
    end
end
s = rank(full(W));
deficiency = n-s-l
pause
% gamma0 = sum(Y*diag(D));
% gamma_vec = gamma0*(0.2:0.2:5);
% C_vecs = [];
% for gamma = 0.2:0.2:5
%     x0 = gamma*ones(n,1);
%     [iter,v_vecs,lmda_vecs,mass_infeas,mass_action_infeas,mu_vecs]=...
%         FpIterationPDCO(Y,Ak,1e-6,500,x0);
%     C_vecs  = [C_vecs exp(lmda_vecs(1:m,end))];
% end

figure(1)
imagesc(sign(Ak)); drawnow;

% figure(2)
% title('Equilibrium Concentrations vs \gamma_{0}');
% ylabel('Equilibrium Concentrations of species');
% xlabel('\gamma_{0}');
% for i = 1:m
%     subplot(3,3,i), plot(gamma_vec, C_vecs(i,:));
% end

% C = exp(lmda_vecs);
% title('Concentration paths of species vs iteration of fixed point algorithm');
% ylabel('Concentration of species');
% xlabel('Iteration');
% figure(3)
% for i = 1:m
%     subplot(3,3,i), plot(C(i,:));
% end

C_soln=[];
for t = 1:50
    x_raw = rand(n,1);
    x0 = 0.001 + 0.999*x_raw./(mass'*Y*D*x_raw);
    [iter,v_vecs,lmda_vecs,mass_infeas,mass_action_infeas,mu_vecs]=...
        FpIterationPDCO(Y,Ak,1e-7,1000,x0);
    C_soln  = [C_soln exp(lmda_vecs(1:m,end))];

    title('Equilibrium points for different starting points');
    ylabel('Concentration of species');
    xlabel('Problem number');
    figure(2)
    for i = 1:m
        subplot(5,5,i), plot(C_soln(i,:)); drawnow;
    end
end